!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
recursive subroutine block_generate (current, var_regrid_flag)
!==================================================================================================================================|

   use parameters
   use block_define
   use interface_hilbert_space_filling_2d

   implicit none

   integer(4), intent(in) :: var_regrid_flag
   type(block), pointer, intent(inout) :: current

   integer(4) :: i, i_p, i_block, i_block_n, j, j_p, j_block, j_block_n, start_x, start_y
   integer(4) :: i_neigh, j_neigh, i_neigh_p, j_neigh_p, i_neigh_c, j_neigh_c
!   real(8), dimension(nx_block, ny_block) :: vx_c, vy_c, vz_c, pr_c, te_c, vx_p, vy_p, vz_p, pr_p, te_p
   type(block), pointer :: child, neigh, next, framework_next, pointer_tmp

!----------------------------------------------------------------------------------------------------------------------------------|
!  before regrid to finer mesh, check the level of neighbor of current block.
!  regrid the neighbor block if the level of neighbor is lower the current block.
!----------------------------------------------------------------------------------------------------------------------------------|
   do j_neigh = -1, 1
   do i_neigh = -1, 1
      neigh => current%neigh(i_neigh, j_neigh)%pblock
      if (associated (neigh)) then
         if (neigh%lv .lt. current%lv .and. neigh%lv .ge. 1) then
            call block_generate (neigh, neigh%cpu_id)
            next => neigh%next
            framework_next => neigh%framework_next
            pointer_tmp => neigh
            call hilbert_space_filling_2d (neigh, pointer_tmp, neigh%start_pos, neigh%end_pos, neigh%start_point,                  &
               ratio_order, neigh%cpu_id)
            if (neigh%cpu_id .eq. cpu_id) then
               next%last => pointer_tmp
               pointer_tmp%next => next
            endif
            framework_next%framework_last => pointer_tmp
            pointer_tmp%framework_next => framework_next
         endif
      endif
   enddo
   enddo

   allocate (current%child(ratio, ratio))

   current%nc_x = ratio
   current%nc_y = ratio
   current%nc = ratio * ratio

   do j_block = 1, ratio
   do i_block = 1, ratio
      allocate (current%child(i_block, j_block)%pblock)
   enddo
   enddo

   if (var_regrid_flag .eq. cpu_id) n_blocks(cpu_id_p1) = n_blocks(cpu_id_p1) + current%nc

!----------------------------------------------------------------------------------------------------------------------------------|
!  interpolation
!----------------------------------------------------------------------------------------------------------------------------------|
   do j_block = 1, ratio
   do i_block = 1, ratio
      child => current%child(i_block, j_block)%pblock

      allocate (child%x(nx_block), child%y(ny_block),                                                                              &
         child%p_pos(n_dims), child%start_pos(n_dims), child%end_pos(n_dims),                                                      &
         child%start_point(n_dims), child%id(n_dims, n_levels))

      child%lv = current%lv + 1
      child%cpu_id = var_regrid_flag
      child%p_pos(1) = i_block
      child%p_pos(2) = j_block
      child%id = 0
      child%id(:, 1 : current%lv) = current%id
      child%nx = nx_block
      child%ny = ny_block
      child%id(:, child%lv) = child%p_pos
      child%dx = current%dx / ratio
      child%dy = current%dy / ratio
      child%dt = current%dt / ratio
      child%t = current%t
      child%parent => current

      start_x = (ng + 1) * (ratio - 1) + (child%p_pos(1) - 1) * nx_no_bnd
      do i = 1, nx_block
         i_p = (start_x + i) / ratio
         child%x(i) = current%x(i_p) + (mod (start_x + i, ratio) - (ratio - 1) / 2.0d0) * child%dx 
      enddo

      start_y = (ng + 1) * (ratio - 1) + (child%p_pos(2) - 1) * ny_no_bnd
      do j = 1, ny_block
         j_p = (start_y + j) / ratio
         child%y(j) = current%y(j_p) + (mod (start_y + j, ratio) - (ratio - 1) / 2.0d0) * child%dy
      enddo

      if (var_regrid_flag .eq. cpu_id) then
         child%var_allocate_flag = .true.

         allocate(child%ro(nx_block, ny_block), child%mx(nx_block, ny_block), child%my(nx_block, ny_block),                        &
                  child%mz(nx_block, ny_block), child%bx(nx_block, ny_block), child%by(nx_block, ny_block),                        &
                  child%bz(nx_block, ny_block), child%en(nx_block, ny_block), child%gx(nx_block, ny_block),                        &
                  child%gy(nx_block, ny_block), child%po(nx_block, ny_block), child%az(nx_block, ny_block))
         allocate(child%fxro_bnd(ny_block, n_neighs), child%fxmx_bnd(ny_block, n_neighs),                                          &
                  child%fxmy_bnd(ny_block, n_neighs), child%fxmz_bnd(ny_block, n_neighs),                                          &
                  child%fxbx_bnd(ny_block, n_neighs), child%fxby_bnd(ny_block, n_neighs),                                          &
                  child%fxbz_bnd(ny_block, n_neighs), child%fxen_bnd(ny_block, n_neighs),                                          &
                  child%fyro_bnd(nx_block, n_neighs), child%fymx_bnd(nx_block, n_neighs),                                          &
                  child%fymy_bnd(nx_block, n_neighs), child%fymz_bnd(nx_block, n_neighs),                                          &
                  child%fybx_bnd(nx_block, n_neighs), child%fyby_bnd(nx_block, n_neighs),                                          &
                  child%fybz_bnd(nx_block, n_neighs), child%fyen_bnd(nx_block, n_neighs),                                          &
                  child%fxpo_bnd(ny_block, n_neighs), child%fypo_bnd(nx_block, n_neighs))

         start_x = (ng + 1) * (ratio - 1) + (child%p_pos(1) - 1) * nx_no_bnd
         start_y = (ng + 1) * (ratio - 1) + (child%p_pos(2) - 1) * ny_no_bnd

         if (ns .eq. 1 .and. interpolation_flag .eq. 1) then
            call model_variables(child%ro, child%mx, child%my, child%mz, child%bx, child%by, child%bz, child%en,                   &
               child%gx, child%gy, child%po, child%x, child%y, nx_block, ny_block)
            call block_interpolate(child%az, nx_block, ny_block, current%az, nx_block, ny_block, start_x, start_y, ratio)
            if (ini_value_save_flag .ne. 0) then
               allocate (child%ro0(nx_block, ny_block), child%mx0(nx_block, ny_block), child%my0(nx_block, ny_block),              &
                  child%mz0(nx_block, ny_block), child%bx0(nx_block, ny_block), child%by0(nx_block, ny_block),                     &
                  child%bz0(nx_block, ny_block), child%en0(nx_block, ny_block), child%az0(nx_block, ny_block))
               do j = 1, ny_block
               do i = 1, nx_block
                  child%ro0(i, j) = child%ro(i, j)
                  child%mx0(i, j) = child%mx(i, j)
                  child%my0(i, j) = child%my(i, j)
                  child%mz0(i, j) = child%mz(i, j)
                  child%bx0(i, j) = child%bx(i, j)
                  child%by0(i, j) = child%by(i, j)
                  child%bz0(i, j) = child%bz(i, j)
                  child%en0(i, j) = child%en(i, j)
                  child%az0(i, j) = child%az(i, j)
               enddo
               enddo
            endif
         else
!            call variables_transform (vx_p, vy_p, vz_p, pr_p, te_p,                                                                &
!               current%ro, current%mx, current%my, current%mz, current%bx, current%by, current%bz, current%en,                     &
!               gm, nx_block, ny_block, 1)
            call block_interpolate (child%ro, nx_block, ny_block, current%ro, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%mx, nx_block, ny_block, current%mx, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%my, nx_block, ny_block, current%my, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%mz, nx_block, ny_block, current%mz, nx_block, ny_block, start_x, start_y, ratio)
!            call block_interpolate (vx_c, nx_block, ny_block, vx_p, nx_block, ny_block, start_x, start_y, ratio)
!            call block_interpolate (vy_c, nx_block, ny_block, vy_p, nx_block, ny_block, start_x, start_y, ratio)
!            call block_interpolate (vz_c, nx_block, ny_block, vz_p, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%bx, nx_block, ny_block, current%bx, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%by, nx_block, ny_block, current%by, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%bz, nx_block, ny_block, current%bz, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%en, nx_block, ny_block, current%en, nx_block, ny_block, start_x, start_y, ratio)
!            call block_interpolate (te_c, nx_block, ny_block, te_p, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%gx, nx_block, ny_block, current%gx, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%gy, nx_block, ny_block, current%gy, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%po, nx_block, ny_block, current%po, nx_block, ny_block, start_x, start_y, ratio)
            call block_interpolate (child%az, nx_block, ny_block, current%az, nx_block, ny_block, start_x, start_y, ratio)
!            call variables_transform (vx_c, vy_c, vz_c, pr_c, te_c,                                                                &
!               child%ro, child%mx, child%my, child%mz, child%bx, child%by, child%bz, child%en,                                     &
!               gm, nx_block, ny_block, 2)
            if (ini_value_save_flag .ne. 0) then
               allocate (child%ro0(nx_block, ny_block), child%mx0(nx_block, ny_block), child%my0(nx_block, ny_block),              &
                  child%mz0(nx_block, ny_block), child%bx0(nx_block, ny_block), child%by0(nx_block, ny_block),                     &
                  child%bz0(nx_block, ny_block), child%en0(nx_block, ny_block), child%az0(nx_block, ny_block))
!               call variables_transform (vx_p, vy_p, vz_p, pr_p, te_p,                                                             &
!                  current%ro0, current%mx0, current%my0, current%mz0, current%bx0, current%by0, current%bz0, current%en0,          &
!                  gm, nx_block, ny_block, 1)
               call block_interpolate (child%ro0, nx_block, ny_block, current%ro0, nx_block, ny_block, start_x, start_y, ratio)
               call block_interpolate (child%mx0, nx_block, ny_block, current%mx0, nx_block, ny_block, start_x, start_y, ratio)
               call block_interpolate (child%my0, nx_block, ny_block, current%my0, nx_block, ny_block, start_x, start_y, ratio)
               call block_interpolate (child%mz0, nx_block, ny_block, current%mz0, nx_block, ny_block, start_x, start_y, ratio)
!               call block_interpolate (vx_c, nx_block, ny_block, vx_p, nx_block, ny_block, start_x, start_y, ratio)
!               call block_interpolate (vy_c, nx_block, ny_block, vy_p, nx_block, ny_block, start_x, start_y, ratio)
!               call block_interpolate (vz_c, nx_block, ny_block, vz_p, nx_block, ny_block, start_x, start_y, ratio)
               call block_interpolate (child%bx0, nx_block, ny_block, current%bx0, nx_block, ny_block, start_x, start_y, ratio)
               call block_interpolate (child%by0, nx_block, ny_block, current%by0, nx_block, ny_block, start_x, start_y, ratio)
               call block_interpolate (child%bz0, nx_block, ny_block, current%bz0, nx_block, ny_block, start_x, start_y, ratio)
               call block_interpolate (child%en0, nx_block, ny_block, current%en0, nx_block, ny_block, start_x, start_y, ratio)
!               call block_interpolate (te_c, nx_block, ny_block, te_p, nx_block, ny_block, start_x, start_y, ratio)
               call block_interpolate (child%az0, nx_block, ny_block, current%az0, nx_block, ny_block, start_x, start_y, ratio)
!               call variables_transform (vx_c, vy_c, vz_c, pr_c, te_c,                                                             &
!                  child%ro0, child%mx0, child%my0, child%mz0, child%bx0, child%by0, child%bz0, child%en0,                          &
!                  gm, nx_block, ny_block, 2)
            endif
         endif
      endif

!----------------------------------------------------------------------------------------------------------------------------------|
!     neighbors
!----------------------------------------------------------------------------------------------------------------------------------|
      do j_neigh = -1, 1
         j_block_n = j_block + j_neigh
         do i_neigh = -1, 1
            i_block_n = i_block + i_neigh
            child%neigh(i_neigh, j_neigh)%pblock => null()
            if (i_block_n .ge. 1 .and. i_block_n .le. ratio .and. j_block_n .ge. 1 .and. j_block_n .le. ratio) then
               child%neigh(i_neigh, j_neigh)%pblock => current%child(i_block_n, j_block_n)%pblock
            else
               i_neigh_p = floor(((i_block - 0.5d0) + i_neigh) / ratio)      !     find the parent block of child's neighbors
               j_neigh_p = floor(((j_block - 0.5d0) + j_neigh) / ratio)
               neigh => current%neigh(i_neigh_p, j_neigh_p)%pblock
               if (associated(neigh)) then
                  if (neigh%nc .ne. 0) then
                     i_neigh_c = i_block + i_neigh
                     j_neigh_c = j_block + j_neigh
                     if (i_neigh_c .eq. 0) i_neigh_c = ratio
                     if (j_neigh_c .eq. 0) j_neigh_c = ratio
                     if (i_neigh_c .eq. ratio + 1) i_neigh_c = 1
                     if (j_neigh_c .eq. ratio + 1) j_neigh_c = 1
                     child%neigh(i_neigh, j_neigh)%pblock => neigh%child(i_neigh_c, j_neigh_c)%pblock
                     child%neigh(i_neigh, j_neigh)%pblock%neigh(-i_neigh, -j_neigh)%pblock => child
                  else
                     child%neigh(i_neigh, j_neigh)%pblock => neigh
                  endif
               else
                  child%neigh(i_neigh, j_neigh)%pblock => null()    !  for physical boundary
               endif
            endif
         enddo
      enddo

   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine block_generate
